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ABSTRACT 


Electrical distribution systems are the last stage of electricity delivery, bridging the gap 
between transmission networks and the customer. Compared to transmission systems, 
there has been comparatively little analysis performed on distribution infrastructure. We 
first establish a minimum load shed optimal switching model. The model seeks to ensure 
all electricity demand is met for a given distribution system by minimizing the amount of 
load not met. This is accomplished by closing and opening switches within the network, 
subject to various electrical physics constraints. We initially test our model against a simple 
toy network to verify proper operation before expanding to increasingly complex networks. 
We culminate with testing against a standardized test circuit. 

The model runs in an open-source programming language on a consumer-grade 
lap-top, reading in data from common formats with solvers readily familiar to the 
optimization community. 
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Executive Summary 


Of the 16 critical infrastructure sectors, electric power is identified as being significantly 
important because it enables the operation of the other sectors. Within the realm of electric 
power, electrical distribution systems consist of the substations, transformers, and lines that 
bridge the gap between transmission systems, which originate at power plants such as coal 
or nuclear, and the end-user. Distribution systems are most readily characterized as the 
power poles and lines running between homes, commercial businesses, and neighboring 
municipalities. 

Working in concert with our associates at Los Alamos National Labs (LANL) and under 
the Joint Warfare Analysis Center (JWAC) - NPS research program, we build a linearized 
optimization model in the Python programming language, implementing the mathematical 
program in the Pyomo environment. The original formulation was provided by our LANL 
associates and expanded upon during model development. At its core, the model delivers 
real and reactive power from an originating substation(s) to loads throughout the network, 
adhering to various electrical engineering physics constraints. 

These constraints include flow balance, balance between the power transmission for phases 
within a single line, adhering to the thermal limits of lines, line losses due to resistance and 
reactance, maintaining nominal system operating voltage, and various topology constraints. 
The entire model is driven by an objective function seeking to minimize the load shed 
throughout the network. 

Testing and analysis is initially performed on a test circuit provided by our associates at 
LANL. The circuit is simple enough that basic operation can be checked by hand calculation 
to verify model behavior. We then specifically examine behavior of the topology constraints 
against another test circuit we designed. We observe the model performing proper operation 
of switches within the network to preclude cycle formation while still minimizing load shed. 

Our testing culminates with application of an Institute of Electrical and Electronics Engi¬ 
neers (IEEE) standardized test circuit. The test circuit possesses more robust and dynamic 
electrical characteristics, specifically imbalanced loading and phase variations. Zero load 
shed is expected with the test circuit and we observe none with our model. The downstream 
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goal of our work here is to perform analysis on real-world distribution systems when they 
are subject to various disruptions. As an intermediate step toward that goal we duplicate 
the IEEE network into “a” and “b” sides, connecting them with a single line containing a 
switch. We proceed to remove the substation from the “b” side and observe the resulting 
behavior as the remaining “a” substation attempts to power the entire system. 

We conclude with recommendations for continued work and analysis. Specifically, our 
model can be folded into attacker-defender analyses to inform researchers on vulnerabilities 
within a given distribution system. Furthermore, testing of real-world systems by triggering 
systematic failures of given components can yield insights to distribution system behavior 
when under duress. Our work is unique due to its focus on distribution systems, its adherence 
to electrical engineering mathematics, and results with guaranteed optimality. 
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CHAPTER 1: 
Introduction 


The U.S. government defines critical infrastructure as “systems and assets, whether physical 
or virtual, so vital to the United States that the incapacity or destruction of such systems and 
assets would have a debilitating impact on security, national economic security, national 
public health or safety, or any combination of those matters” (White House, 2013). Critical 
infrastructure systems are those which deliver the energy to heat and cool our homes, the 
water we drink, the communication systems we enjoy, and the transportation systems that are 
integral to the American way of life. Furthermore, these systems serve as the foundation for 
the economy, security, and health of our nation (Department of Homeland Security, 2017). 
Due to the significance of these systems, there is considerable interest in understanding, 
modeling, and predicting of infrastructure system response to disruption. It is of paramount 
interest at the presidential level that we “maintain secure, functioning, and resilient critical 
infrastructure” (White House, 2013). Electric power, one of these critical infrastructure 
systems, is the focus of this research. 

1.1 The Importance of Electric Power 

As defined by the Department of Homeland Security (DHS), there are 16 critical infras¬ 
tructure sectors (Department of Homeland Security, 2017). Of these, electric power is 
particularly important because it provides an “enabling function” for the others in funda¬ 
mental ways (White House, 2013). 

The overall structure of electrical power infrastructure can be loosely categorized into 
three main parts: generation, transmission, and distribution. Generation systems, such as 
solar, hydro, wind, nuclear, or gas turbine power plants, produce electricity. Transmission 
systems move electricity in bulk from the originating power plant to substations and are 
most commonly recognized as the massive towers carrying heavy lines across rural stretches 
of highway. The distribution system exists between the substation and end customer, most 
readily characterized as the telephone poles abundant in any residential and commercial 
area. 
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Disruptions to the electric power system can be caused by a variety of sources. Foremost 
among these are natural disasters such as hurricanes or earthquakes. For example, Hurricane 
Sandy in 2012 caused power losses to over 8,000,000 homes across 17 states (Webley, 2012). 
Other potential causes of electric power disruption are operator error (e.g., the Three Mile 
Island incident; see Nuclear Regulatory Commission, 2013) or the actions of an adversary 
(e.g., the coordinated attacks on the Metcalf substation in 2013 near San Jose, CA; see 
Smith, 2014). Consequences for these disruptions vary depending on the scope and scale of 
the damage as well as the types of electricity infrastructure damaged. A car hitting a power 
pole may interrupt power to a single home, but a worst-case solar event electromagnetic 
pulse (EMP) would take an estimated four to ten years to repair and leave most of the United 
States without power (Federal Energy Regulatory Commission, 2010). 

Regardless of the cause, the consequences of a disruption have significant negative impli¬ 
cations for our welfare, both in terms of loss of life as well as negative economic impacts 
(Anderson and Bell, 2012; Rose et al., 1997). Generation facilities are frequently taken 
online and offline to balance power supply and load demand, but unexpected failure of a 
generation facility can cause problems (Knaus, 2017). Disruptions to transmission systems, 
either from the loss of a high-voltage line or a substation, can impede the ability to move 
bulk power and cause long-distance power outages (Cava, 2017). A disruption to a distri¬ 
bution system causes downstream customers to potentially go without electricity, but these 
effects are often local in nature. 

There is a tactical need to understand the strengths and weaknesses of our own electric grid 
as well as those of potential adversaries. Domestically, identifying weaknesses in our grid 
informs our defense efforts so we can minimize the impact of adversarial action against 
us, because as reporters, politicians, and utility operators know, disruptions are devastating 
(Williams and Bennett, 2016). From an offensive perspective, we desire to know the effects 
of targeting an enemy’s infrastructure to achieve our tactical objective, but we desire to 
restrict the consequence to only what’s necessary to achieve victory. Militarily, the second 
order effects of taking out an entire transmission system in an enemy territory might prove 
too costly. Research has shown that in large scale power outages, ordinary civilians suffer 
and death rates increase (Anderson and Bell, 2012). If we seek to confine power disruption 
to as small an area as possible, attacking the transmission system is off the table. Hart et al. 
(2014) provide a discussion of the tensions and tradeoffs for targeting infrastructure in a 
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military context. 


Prior research and analysis efforts at the Naval Postgraduate School (NPS) have focused 
primarily on the behavior and response of transmission systems. There are multiple re¬ 
search papers and theses on the subject of transmission system modeling (e.g., Rose, 2007; 
Salmeron et al., 2012). Comparatively, little effort has been applied to distribution systems. 

The focus of this thesis is understanding how a distribution system will respond to the loss 
of components (e.g., substations), and in particular the presence and size of outages. In 
response to component failure, a distribution system will adjust its connectivity by opening 
and closing switches to continue powering as much of the grid as possible. We capture this 
behavior in a prescriptive decision Alternating Current Power Flow (ACPF) model which 
changes flow by changing connectivity via switch position adjustments to minimize the 
amount of load shed. 

We implement the ACPF model in the Python programming language (Python Software 
Foundation, 2017), primarily relying on the Pyomo optimization package (ffart et al., 2012, 
2011). We verify model operation by application to a toy model and then expand to a real 
data set from a large Northeastern utility provider. 

1.2 Thesis Structure 

Chapter 2 provides cursory overviews of the academic disciplines underpinning this work. 
Additionally, we review the relevant literature and research pertaining to our analysis. In 
Chapter 3 we introduce the ACPF model and demonstrate its application with a small series 
of toy data sets. Chapter 4 describes our application of the model to a standardized test 
circuit. We conclude in Chapter 5 by summarizing our findings along with recommendations 
for future research. 
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CHAPTER 2: 
Background 


This chapter provides background information on the disciplines involved in this thesis as 
well as pertinent previous research. 

2.1 Technical Background Overview 

Our work encapsulates two major academic domains: mathematical optimization and 
Electrical Engineering (EE). While it is beyond the scope of this thesis to fully inform 
the reader on both these subjects, we seek to provide a brief overview to aid reading our 
work. 

2.1.1 Mathematical Programming 

Mathematical programming involves theory and methods for solving for extreme points, 
either maxima or a minima, of functions on sets defined by linear and/or non-linear con¬ 
straints. The primary objective function for which a maxima or a minima is being found, 
commonly involves a maximum profit or lowest cost. The constraints are linear equality 
and inequalities representing limitations imposed on the system, such as how much money 
can be spent or the capacity limits of a path. For a thorough introduction to optimization 
and mathematical programming, see the textbook Optimization in Operations Research 
(Rardin, 2016). 

2.1.2 Electrical Engineering 

For a deeper understanding of EE as it applies to electric power systems we recommend the 
work Elements of Power Systems (Sadhu and Das, 2015). 

Electrical transmission and distribution systems are frequently 3-phase Alternating Current 
(AC) systems. The three phases are conventionally a, b, or c and a given load or line in the 
network can be any combination of the three phases. The supply and demand quantities 
in the network are measured in watts, with transmission systems generally operating in 
Megawatts (Mw) while distribution systems frequently are measured in Kilowatts (Kw). 
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Generated electricity consists of real and reactive components. Real power produces kinetic 
work (e.g., motors spinning), while reactive power is the electromagnetic hindrance of 
electrical processes. Correspondingly, loads in the network consume real and/or reactive 
power. Individual power lines possess both resistance and reactance, qualities that hinder 
real and reactive power flow, respectively. 


2.2 Previous Research 

2.2.1 Attacker-Defender and Defender-Attacker-Defender Analysis 

Attacker-Defender (AD) analysis is a two-part optimization problem. The attacker seeks to 
maximize the minimum value the defender can obtain from the system post attack. There 
are three key assumptions in AD modeling, “(1) the attacker’s and defender’s actions are 
sequential, (2) the attacker has a perfect model of how the defender will (or should) optimally 
operate the system, even after an attack, and (3) the attacker will manipulate that system to 
his best advantage” (Brown et al., 2006). 

Defender is a bit of a misnomer in an AD model. Nothing is being hardened or made 
invulnerable to attack; in reality the defender is simply a system user who best operates post¬ 
attack. However, it is possible to expand a two-level AD model into a three-level Defender- 
Attacker-Defender (DAD) model, where the system is initially hardened to preclude or 
mitigate the most damaging attack (Brown et al., 2006). Afterward and with knowledge of 
the defensive action, the attacker performs a most damaging attack, and lastly the system 
operator works to obtain best system performance. 

Salmeron et al. (2004) applies attacker-defender optimization supported by a heuristic to 
an electric power model. They account for a great deal of electrical physics considerations, 
abstaining from the simplifications of a pure flow model. Their results demonstrate the 
ability to identify the key components in a network that, if attacked, would be most crippling 
to operation and simultaneously represent the greatest need for hardening and resiliency 
efforts. Analysis is performed on Institute of Electrical and Electronics Engineers (IEEE) 
Reliability Test System (RTS) 96. 

Salmeron et al. (2005) introduce a Graphical User Interface (GUI) for electric power 
infrastructure analysis, named Vulnerability of Electric Grids Analyzer (VEGA), which 
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is supported by an optimal power flow linear program and an interdiction model based 
on heuristic algorithms. The heuristic algorithm is used to “find the combination of 
facilities (generating units, buses, transmission lines, transformers and/or substations) whose 
destruction would cause the largest or nearly largest disruption in the system at a given 
point in time, which in turn suggests the system components most in need of hardening” 
(Salmeron et al., 2005). Tests and analysis are performed on two real world transmission 
systems. Additionally, a time-series component is incorporated with acknowledgement of 
peak loading time-frames being the most devastating windows of attack. 

Salmeron and Wood (2009) extend prior work to formulate and solve DAD problems for 
electrical transmission grids, identifying critical components sensitive to terrorist attack. 
Applying a generalized version of Benders decomposition, Salmeron et al. (2009) identify a 
set of components whose interdiction causes maximum economic loss to customers. Their 
test problems “describe a regional power grid in the United States with approximately 5000 
buses, 6000 lines, and 500 generators” (Salmeron et al., 2009). 

Salmeron et al. (2011) perform DAD analysis on Vandenberg Air Force base and its sur¬ 
rounding area. Their efforts identify the vulnerable components whose destruction would 
be maximally disruptive to system operation and presume adversaries are aware of those 
defensive preparations, in turn making maximally destructive attacks. Lastly, the system 
is optimally operated post-attack. This focus on a Department of Defense (DOD) entity is 
unique because the surrounding area benefits from excess generation capacity. 

Further modeling efforts focus on the Guam Power Authority (GPA), identifying critical 
sets of components whose destruction would most impede GPA’s ability to provide power 
(Salmeron et al., 2012). Subsequently, these sets of components benefit the most from 
hardening efforts. 

2.2.2 Previous NPS Theses 

Analysis on electric power systems using AD and DAD models has been driven and sup¬ 
ported by many student theses. 

Stathakos (2003) performs research contributing to the development of VEGA. Alvarez 
(2004) demonstrates empirical speed ups to Benders decomposition on electric grid inter- 
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diction problems, tested against IEEE RTS networks. Optimality gap reductions occurred 
alongside speed improvements. 

Carnal (2005) further improves on VEGA by inclusion of the Xpres-MP software package, 
reducing run-time by 65% - 85%. 

Schneider (2005) performs DAD analysis on the eastern interconnection transmission grid 
and determine hardening levels necessary to reduce risk below acceptable levels. Jonart 
(2008) does follow-on work, making algorithmic improvements, yielding results within 1% 
of optimality. 

Ang (2006) performs a Mixed-integer program (MIP) optimization of electric grid recovery 
efforts after natural disaster or terrorist attack. Constraints include availability of repair 
resources, penalties for unserved electricity demand, and repair time planning horizons. 
Testing is performed on an IEEE 300 bus transmission network. 

Rose (2007) develops a mathematical model to identify critical components to defend in an 
electric grid, given limited resources. Using a pure-flow model for the electric grid. Rose 
reduces solution times from previous efforts while simultaneously finding solutions that are 
40% better. 

2.2.3 Damage Anticipation Modeling on Electrical Systems 

Guikema et al. (2010) use regression and data mining techniques to estimate required 
repair time and replacement utility poles and transformers in the wake of electric power 
infrastructure damage resulting from natural disaster, notably hurricanes. Their efforts 
are predictive in nature and seek to best prepare responders for repairs after a disruption. 
Ouyang and Duenas-Osorio (2014) build a probabilistic model to determine economic, 
social, technical, and organizational resilience to hurricane damage. 

Maliszewski and Perrings (2012) use regression modeling to build spatial outage visual¬ 
izations and identify geographic factors with strong correlation to minimal outage risks in 
Phoenix’s electrical distribution system. They find a strong correlation of neighborhoods 
near high priority customers such as hospitals experiencing shortened outages because of 
“the quality of infrastructure provided and the priority given to restoring supply to those 
users” (Maliszewski and Perrings, 2012). Proximity to vegetation was statistically signifi- 
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cant as well while other services, such as police and fire stations or schools, all which were 
insignificant or negatively correlated. 

A variety of natural disasters, notably the 2010 earthquakes in Haiti and Chile, the 2011 
earthquake and tsunami in Japan, as well as multiple hurricanes, have served as the sources 
of data for various resiliency modeling efforts. MacKenzie and Barker (2012) build a re¬ 
gression model and mixed effects model to anticipate rates of recovery post-disruption and 
anticipate a system’s resiliency. Using data from Hurricane Ike, Wei et al. (2013) develop a 
“spatial-temporal non-stationary random process to model large-scale disruptions of power 
distribution induced by severe weather.” Their stochastic approach yields resiliency esti¬ 
mates for components and identifies dynamic neighborhoods within the network. Kwasinski 
(2016) performs data analysis and regression modeling to inform future damage anticipation 
efforts and design considerations for distribution systems. 

2.2.4 Optimal Switching and RTS usage 

Delgadillo et al. (2010) construct an optimization approach to defend a grid against terrorist 
attack by calculating the best operator response (e.g., switching) to disruption. Their 
work is a two-part optimization, first identifying the most damaging terrorist attack, 
measured in Mw not delivered, and then identify the optimal switch realignment in 
response to that attack. 

Additionally, de Assis et al. (2015) consider optimization of switch allocation for distribution 
systems. They use a distribution system and identify the ideal locations for switches under 
various system design constraints. Specifically, their research minimizes the cost of switch 
allocation, a system design consideration. They also minimize energy not supplied while 
meeting various regulatory agency reliability requirements and system capacity necessities. 
Their solutions are determined by “memetic algorithms, supported by well-established 
genetic algorithm concepts and local search optimization procedures.” 

The IEEE 118-bus and the RTS 96 test system are studied for optimal transmission system 
switching with a mixed integer program (Hedman et al., 2009). Their model demonstrates 
the test networks can be operated to meet N-l secure standards while simultaneously 
cutting costs. Another instance of RTS usage is the work of Alguacil et al. (2003) in their 
transmission expansion planning optimization problem. 
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CHAPTER 3: 
Model 


We develop an alternating current (AC) model for electrical distribution system that in¬ 
tegrates switching and adherence to full 3-phase AC mathematics. Instead of suggesting 
likelihoods of what could happen in a disaster or terrorist scenario, we model what will 
happen for any disruption or series of disruptions. 

3.1 Model Formulation 

Our Alternating Current Power Flow (ACPF) model is built in collaboration with our 
partners at Los Alamos National Laboratory (LANL), notably Harsha Nagarajan and Arthur 
Barnes, with significant sourcing from their previous work (e.g., Barnes et al., 2017; 
Nagarajan et al., 2016). 

The entirety of our ACPF model is presented in NPS format for a given Power Grid (PG). 
Discussion of the objective function and constraints follows in section 3.1.2. 

3.1.1 Mathematical Formulation 

Sets and Indices 


i g N 
ScNxN 
NT c <5 
SQ NT 
3K c NxN 

P - {a, b, c} 


set of nodes (buses) in PG (alias j) 

set of undirected edges (lines and transformers) in PG 

set of edges without transformers in PG 

set of edges with switches in PG 

set of directed arcs (lines and transformers) in PG 

(l j) e £ (6 j) 6 £ A (j, i) G £ 

set of phases in PG 
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k e Pi c <P set of phases on bus i which contains a load or a generator 


Pij c P set of phases on line (ij) in PG 
C set of cycles in PG 
N(C) set of nodes in cycle C 
<5(C) set of edges in cycle C 

Parameters [units] 

Api k cost for shedding real power at bus i on phase k [$/watt] 

Aqi k cost for shedding reactive power at bus i on phase k [$/watt] 

Tjj£ thermal limit on edge (i.j) for phase k [Amps] 
vj ,, v u , voltage lower and upper bounds at bus i for phase k [Kv] 
i imaginary number constant 

dpj k + i dcp k AC power demand at bus i on phase k [Kw and Kvar] 

gp u ik + i gq 1 ' k existing AC power generation capacity at bus i on phase k [Kw and Kvar] 

Rjj'k + i Xjj,k impedance of edge ( i,j ) on phase k in PG [Q] 

Zij impedance of edge (IJ). defined as: 

Z,J = z,je^-, Z u = Zije-W where, Z,j = R,j +iX l j. Z,, = R„+iX tJ [01 
M u Mij = (1.05V” OW ) 2 - (0.95P""'") 2 V/./ 6 6 [Volts 2 ] 

Pg allowed variation in flow between phases on edge (ij) [0.15 for transformer 
edges, 1.0 for all other edges] 

yn° m nom j na ] voltage the system is operating at [Kv] 
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Decision Variables [units, when appropriate] 


vik voltage at bus i on phase k [Kv] 

Xij determines if directed arc (i,j) carries power flow in PG [binary] 
Pij k real AC power flow on arc (i, j ) [Kw] 
qj h k reactive AC power flow on arc (i,j) [Kvar] 

Ipi'U + i lqi,k AC power load that is shed at bus i on phase k 


gpi^k + i gqi,k existing AC power generation capacity at bus i on phase k 


Formulation 

min I (Api,klPi,k + Mi,klquk ) (1) 

i£N,k£Pj 

S.t. 


Power flow constraints 


9Pi,k - dpi,k + lpi,k = J] Pij,k - J] PiUk 

ijeS ji&£ 

gqik — dqi k + lqi,k — ^ ^ qij,k ~ ^ ] qji.k 

ij£& ji£& 


(l-Av) 


YjkePij Pij,k 


< Pij,k ^ (1 + fiij) 


JukePi, Pij,k 


— qij,k — (1 + Ay) 


TjkePtj qij,k 


Vi € ft, k £ Pi 


Vi £ ft, k £ Pj 


Vi j £ ft, k e Pij 


Vi j £ ft, k £ Pij 

Vi j £ ft, k £ Pi 
Vi j £ ft, k £ Pi 


a - Ai) s S(U Ay) v,7 ao Ay 

Pi M + tijjk * *jTfj*<yrf +3 Vi/ 6 ,*,* e *>y 

Pij.* + dij,k < XijTij, k V; iom yJ 6 V/y 6 6 Pij 

( v j,a ) 2 - (Pi,a)~ + 2 (RPpij'C + XPqipa + Pij^ + Xjj qij'b + Rjj Pij,c + Xjj qij,c) 
>-Mij(\-Xij) Vij £ NT 

( Vj,b f - (v ub ) 2 + 2(R^pij b + xjfqij'b + R^j Pij,c + x’i-qij'C + R-j'p^a + X^q ijA ) 

> — Xjj) Vij £ NT 


(2a) 

(2b) 

(3a) 

(3b) 

(4) 

(4a) 

(5a) 

(5b) 
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(5c) 


(vj,cY - (Vi, c y + 2 {R C ijPij, c + x u<lij,c + RjJ Pij,a + X]Jq iM + R--Pij,b + Ku<lij,b) 


rcb, 


> -M,j(\ — Xij) Vij £ NT 

(Vj,a) 2 - (Vi,a) 2 + 2 {R^Pij'd + X^ a <Jij,a + Rjj Pij,b + X ij + Ejf Pij,c + X^j <Jij,c) 
< Mij{\ - Xij) V ij g NT 

( v j,b) 2 ~ (Vi, b y + 2 {R^Pij'b + xffqijj, + Rjj pij, c + Xjj qg ;C + R^pg, a + Kg'<lij,a) 


< Mg{\ -Xij) Vij £ NT 

( Vj,c ) 2 _ ( v i,c) 2 + 2 {RgPij,c + Xfj C qij tC + Rjj Piij >a + Xjj qjj^a + RjjjPij,b + Xdg(pj,b) 


< Mjj{ 1 - Xij) 

(0.95 v; wm f < vf k < (i.05v; iom ) 2 
0 < qp lM < gp u i k 

gquk < gq\[ k 
- (R/i.k < gq“ k 


Vi j £ NT 
Vi g N,k e Vj 
Vi G N, k G Pi 
Vi G N, k G Pi 
Vi G N, k G Pi 


lpi, k < dp kk 
lqi, k — dq lM 


Vi &N,k G Pj 
Vi £ N,k £ Pi 


(5d) 

(5e) 

(5f) 

( 6 ) 

(7a) 

(7b) 

(7c) 

(8a) 

(8b) 


Topology constraints 

y xij<\N{o\-\ 

ijeS(C) 

Xij + Xjj < 1 
Xij + Xjj > 1 


VC (9) 

V(iJ)eS (10a) 

V(i,j)£&\S (10b) 


3.1.2 Model Discussion 

The objective function (1) minimizes the cost associated with all real and reactive load 
shed due to the loss of generation capacity at bus i. The cost parameters, Api, k and Aq lk 
(nominally set to 1) scale the total value of Kw and Kilovolt-ampere reactive (Kvar) shed 
by the model. 

Constraints (2a) and (2b) maintain balance of flow throughout the network, node by node, 
between generation, load, and load shedding, for real and reactive power, respectively. 
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Constraints (3a) and (3b) ensure that real and reactive power flows for each phase are within 
stated tolerances. A large difference in flow on a single line between two phases isn’t 
desirable from an EE perspective. The fiij terms are tolerance parameters: the higher the 
Pij value, the wider the tolerance between phases that is allowed. 

Constraint (4) ensures the thermal limits of line (z, j) are not exceeded. These limits are 
enforced only if a line is being used (i.e., Xjj = 1). Constraint (4a) is a linear version to 
allow each x-,j to control whether flow is permitted on arc (z, j). 

Constraints (5a) through (5f) account for flow losses across all possible phase combinations 
on line (z j) due to resistance and reactance values for that line. Again, these constraints are 
enforced only if a line is being used (i.e., x,y = 1). 

Constraint (6) maintains voltage at node z within 95% to 105% of the nominal voltage for 
that node, as mapped by its location from the nearest upstream transformer. 

Constraints (7a), (7b), and (7c) ensure nodes only generate real and/or reactive power if 
capable, as appropriate. 

Constraints (8a) and (8b) only allow load to be shed at nodes with real or reactive load. 
This prevents the formulation from spuriously shedding load at nodes without load, simply 
to satisfy constraint mathematics. 

Constraint (9) precludes cycles in the network by ensuring the number of active edges in 
a cycle is one less than the number of edges in that cycle. This constraint is implemented 
differently than the others. Initially, we solve the ACPF model without the cycle breaking 
constraint active. Once a solution is obtained we analyze it to identify cycles within that 
solution. For each cycle in the incumbent solution, a copy of constraint (9) is added and 
then the ACPF model is solved again. 

Constraint (10a) ensures there is flow in at most one direction (namely, on one arc) for a 
given edge, and constraint (10b) ensures that edges without a switch are active in at least 
one direction. If the switch on edge (z, j) is open, then both x,j and Xj l will be zero. 
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3.2 Model Implementation 

The model is built in the Python programming language, an interpreted, high-level, general 
purpose programming language (Python Software Foundation, 2017). The modular and 
flexible nature of Python enabled us to integrate all computational tasks within a single pro¬ 
gram. The data read and processing steps use several standard Python packages, including 
j son, numpy, and math packages, as well as networkx (Hagberg et al., 2008) for analysis. 
The mathematical program is implemented in Pyomo, a “Python-based, open-source op¬ 
timization modeling language with a diverse set of optimization capabilities” (Hart et al., 
2012 , 2011 ). 

Our test networks and the real life network are all in JavaScript Object Notation (JSON) 
format (Ecma International, 2017), with three dictionaries in the data file, one for the nodes, 
one for the edges in the network, and another for more detailed information on the nodes. 
JSON, consisting of name and value pairs for its data storage is a natural fit for the Python 
(and Pyomo) objects we use to construct the model. 

After constructing all of the Python storage objects necessary for the model we have two 
primary loops to prepare the data, one for the nodes and one for the edges. During these 
loops we read in all of the requisite details for the distribution system. First, we read in node 
details including attributes such as the type of node, e.g., substation, shunt capacitor, or 
load, a unique node identification number, and specific data items such as dp^k, dq^k, gp u k , 
v l .,, and more. Next, we read in data for the lines (i.e., edges) in the system. We capture the 
Rjj^k and for lines, which nodes in the network a given line connects, and the number 
of phases in a line. For lines that are transformers we capture the voltage difference across 
the transformer and map these voltages through the network as V[ wm . 

We store all data in Python dictionaries. Dictionaries are an ideal storage medium for the 
network data because their structure uses key-value pairs and in this case we use the unique 
node identification number to index node data, and the (tail, head) node identification pairs 
for an edge to index edge data. Furthermore, Python dictionaries naturally integrate with the 
Pyomo optimization package for forming Pyomo objects. For example, once a dictionary 
is created with Rjj,k, it can be fed directly to the Pyomo resistance object. 

Once the data is read in and prepared, we build the ACPF model in Pyomo and solve it with 
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the CPLEX solver (IBM, 2017). 


3.3 Model Testing 

We begin with a basic test circuit that has no cycles or switches, a minor design simplification 
to facilitate incremental model construction. A depiction of the basic test circuit appears in 
Figure 3.1, with associated parameters in Table 3.1 and Table 3.2. 


Substation (3,0) 

t 


Substation 

Transformer 



(TO) (3,5) 


Load (3,139) 



Service 

Transformer 



( 3 , 119 ) 


( 3 , 139 ) 


I 

Shunt Cap Bank (3,139) 


The Basic Test Circuit consists of a substation with functionally infinite capacity, 
a transformer, connective line leading to a service transformer, and a commer¬ 
cial grade customer with moderate loading and a shunt capacitor bank. When 
operated, no load should be shed in this test network. 

Figure 3.1. Basic Test Circuit 


3.4 Model Performance on Basic Test Circuit 

The ACPF model, running on an early 2015 MacBook Pro with a 2.7 GHz i5 processor, 
solves the basic test circuit in less than a second. The objective function (1) returns a result 
of 0.00 corresponding to zero load shed, as anticipated. We observe pij^ flows from the 
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Table 3.1. Basic Test Circuit Node Data 


Node 

Type 

9Pi,k 

c Kk 

dpi,k 

dqi.k 

1 ynom 

(3,0) 

Substation 

OO 

OO 

0 

0 

69 

(3,5) 

Connection 

0 

0 

0 

0 

7.62 

(3,119) 

Connection 

0 

0 

0 

0 

7.62 

(3,139) 

Load, 

Shunt 

capacitor 

0 

400 

57.55 (a) 
55.48 (b) 
60.65 (c) 

29.29 (a) 
27.51 (b) 
30.78 (c) 

0.277 


Load for each phase (a-c) indicated explicitly 


Table 3.2. Basic Test Circuit Edge Data 


Edge 

Type 

Phases 

Tij,k 

R q,k 

Xii,k 

(3,0),(3,5) 

Transformer 

a , b, c 

120.773 

0.0217 

0.3548 

(3,5),(3,119) 

Line 

a, b, c 

498.75 

0.1334 (a) 
0.1354(b) 
0.1343 (c) 

0.3124 (a) 
0.3053 (b) 
0.3089 (c) 

(3,119),(3,139) 

Line 

a, b, c 

32.808 

0.0009 

0.0053 


Rjjj, and Xij.k for each phase listed explicitly when variations exist. 


substation to the load of 57.55, 55.48, and 60.65 Kw for phases a,b, and c, respectively, 
values matching the demands in node (3,139). 

Due to the presence of a shunt capacitor bank on node (3,139) all reactive demand is 
satisfied locally at the commercial customer, resulting in no qij^ flow. Iterative testing, 
namely removal of the shunt capacitor bank, causes qq^ flow (in Kvar) satisfying all dq i k . 
This behavior confirms proper ACPF model operation regarding reactive power flow as 
well. 

At this point our ACPF model has performed as expected on the basic test circuit, success¬ 
fully operating while adhering to power flow constraints (2a) through (8b). 

3.5 Model Testing With Cycles 

Once satisfied with operation on the basic test circuit we added switching controls and the 
cycle breaking constraint (9) to complete implementation of the formulation in our ACPF 
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model. At this point we build a basic test circuit with a cycle and switches for testing, as 
depicted in Figure 3.2, with associated data in Table 3.3 and Table 3.4. 



Switch 2,4 



The Basic Test Circuit With Cycle consists of a substation with functionally infinite 
capacity and a moderate load, each with their own transformer. Between them is 
a cycle of four nodes, with two of the four edges containing switches. Unlike the 
Basic Test Circuit, there is no shunt capacitor at the load. 

Figure 3.2. Basic Test Circuit with Cycle 


Table 3.3. Basic Test Circuit with Cycle Node Data 


Node 

Type 

9Plk 


dpi,k 

d(Ji,k 

ynom 

Sub 

Substation 

oo 

oo 

0 

0 

100 

1 

Connection 

0 

0 

0 

0 

10 

2 

Connection 

0 

0 

0 

0 

10 

3 

Connection 

0 

0 

0 

0 

10 

4 

Connection 

0 

0 

0 

0 

10 

Load 

Load 

0 

0 

30.0 

20.0 

1 


Loading for each phase a - c is identical. 
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Table 3.4. Basic Test Circuit with Cycle Edge Data 


Edge 

Type 

Phases 

Tij,k 

R ij,k 

*ij,k 

Sub,l 

Transformer 

a, b, c 

100 

0.0217 

0.3548 

1,2 

Line 

a, b, c 

100 

0.1334 

0.3124 

1,3 

Line 

a, b, c 

100 

0.1334 

0.3124 

2,4 

Line 

a, b, c 

100 

0.1334 

0.3053 

3,4 

Line 

a, b, c 

100 

0.1343 

0.3089 

4,Load 

Transformer 

a, b, c 

100 

0.0009 

0.0053 


Rjj.k and X/jk for each phase listed explicitly when variations exist. 


The Basic Test Circuit with Cycle contains realistic R/j and Xg values, drawn from real 
networks. Possessing a single load with its own service transformer and a substation with 
functionally infinite gp t ,k and gq t .k, the basic test circuit with cycle is electrically sound and 
designed to shed no load when solved. 


3.6 Model Performance on Basic Test Circuit with Cycle 

Our model solves the basic test circuit with cycle initially, not enforcing the cycle breaking 
constraint, resulting in the behavior summarized in Table 3.5. 


Table 3.5. Basic Test Circuit with Cycle Results before Cycle Breaking 
Constraint 


Edge 

Type 

Switch / Position 

x ij 

PijX 

Qij-k 

Sub,l 

Transformer 

Lalse, N/A 

1 

30.0 

20.0 

1,2 

Line 

Lalse, N/A 

1 

15.0 

10.0 

1,3 

Line 

True, Closed 

1 

15.0 

10.0 

2,4 

Line 

True, Closed 

1 

15.0 

10.0 

3,4 

Line 

Lalse, N/A 

1 

15.0 

10.0 

4,Load 

Transformer 

Lalse, N/A 

1 

30.0 

20.0 


Pij,k and qij t k flows the same for all phases. 
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At this point, we check the incumbent solution for cycles (using the networkx package) 
and identify the cycle consisting of nodes 1, 2, 3, and 4. The cycle breaking constraint (9) 
is enforced and the ACPF model runs again, this time producing a solution with no cycles 
and zero load shed, as summarized in Table 3.6. Cycle breaking is restricted to opening a 
switch on lines that have switches. 


Table 3.6. Basic Test Circuit with Cycle Results after Cycle Breaking Con¬ 
straint 


Edge 

Type 

Switch / Position 

x ij 

Pij,k 

Qij,k 

Sub,l 

Transformer 

False, N/A 

1 

30.0 

20.0 

1.2 

Line 

False, N/A 

1 

0.0 

0.0 

1,3 

Line 

True, Closed 

1 

30.0 

20.0 

2,4 

Line 

True, Open 

0 

0.0 

0.0 

3,4 

Line 

False, N/A 

1 

30.0 

20.0 

4,Load 

Transformer 

False, N/A 

1 

30.0 

20.0 


Pij,k and qtj'k flows the same for all phases. 


As can be observed, the cycle breaking constraint has forced the switch on line 2,4 to open. 
Flow redirected to meet demand while avoiding a cycle being present in the network. 

3.7 Model Validation 

Having tested our ACPF model’s basic functions we test it against an IEEE RTS network for 
validation. The IEEE RTS is a standardized test circuit whose behavior is well understood 
by researchers and academia. Furthermore, phase changes, multiple branching paths, and 
uneven loading make it a more robust test case than the basic test networks used thus far. 
Observing proper behavior of the ACPF model on the RTS network serves as validation to 
our modeling and these efforts are described in Chapter 4. 
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CHAPTER 4: 
Results 


4.1 IEEE RTS 13-bus feeder 

We choose the IEEE RTS 13 node test feeder because it is relatively small and frequently 
used to test analysis software (IEEE Power and Energy Society, 2017). Furthermore, the 
network itself mimics real world distribution and is “characterized by being short, relatively 
highly loaded, a single voltage regulator at the substation, overhead and underground lines, 
shunt capacitors, an in-line transformer, and unbalanced loading” (IEEE Power and Energy 
Society, 2017). A network with these characteristics will test many facets of our model 
in a more robust capacity than the test networks in Chapter 3. The network is constructed 
without the potential for cycles however, so our cycle breaking constraint is not exercised 
when solving the IEEE RTS 13 node test feeder. The circuit is depicted in Figure 4.1. 



Figure 4.1. IEEE 13 Node Test Feeder. The substation is at node 634 
with its transformer between node 634 and 633. A voltage regulator exists 
between node 632 and 650 but has no impact on network behavior within 
our model because node 650 is an unloaded bus. 
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Full electrical characteristics for the 13 node test feeder can be found on the IEEE Power and 
Energy Society website (IEEE Power and Energy Society, 2017). The basic characteristics 
and those specifically pertinent to our ACPF model are described in Tables 4.1 and 4.2. 
Additionally, there is a distributed load in the IEEE 13 node test feeder along the line 
between nodes 632 and 671. We capture this load as an extra node named “dload,” existing 
equidistant between nodes 632 and 671. The distributed load node has dp kk values of 10, 
38, and 68, with dq kk values of 17, 66, and 117, all for phases a, b, and c, respectively. 


Table 4.1. IEEE 13 Node Test Feeder Node Data 


Node 

Type 

gplk 

9<t Lk 

dp L k 

dcp.k 

ynom 

634 

Substation, 

Load 

oo 

oo 

160(a) 

120(b,c) 

110(a) 

90(b,c) 

115 

633 

Connection 

0 

0 

0 

0 

4.16 

632 

Connection 

0 

0 

0 

0 

4.16 

650 

Connection 

0 

0 

0 

0 

4.16 

645 

Load 

0 

0 

170(b) 

125(b) 

4.16 

646 

Load 

0 

0 

230(b) 

132(b) 

4.16 

671 

Load 

0 

0 

385(a,b,c) 

220(a,b,c) 

4.16 

684 

Connection 

0 

0 

0 

0 

4.16 

652 

Load 

0 

0 

128(a) 

86(a) 

4.16 

611 

Load, 

Shunt 

capacitor 

0 

100(c) 

170(c) 

80(c) 

4.16 

680 

Connection 

0 

0 

0 

0 

4.16 

692 

Load 

0 

0 

170(c) 

151(c) 

4.16 

675 

Load, 

Shunt 

capacitor 

0 

200 

(a,b,c) 

485(a) 

68(b),290(c) 

190(a) 

60(b),212(c) 

4.16 


Load for each phase (a-c) indicated explicitly 
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Table 4.2. IEEE 13 Node Test Feeder Edge Data 


Edge 

Type 

Phases 

Tij,k 

634,633 

Transformer 

a, b, c 

1500 

633,632 

Line 

a, b, c 

340 

632,650 

Line 

a, b, c 

730 

632,645 

Line 

b, c 

230 

645,646 

Line 

b, c 

230 

632,dload 

Line 

a, b, c 

730 

dload,671 

Line 

a, b, c 

730 

671,684 

Line 

a, c 

230 

684,611 

Line 

c 

230 

684,652 

Line 

a 

165 

671,680 

Line 

a, b, c 

730 

671,692 

Line 

a, b, c 

1500 

692,675 

Line 

a, b, c 

260 


Rjj j, and Xjj j. for the network are 3 x 3 matrices. We refer our readers to the IEEE RTS 

reference materials for specific values. 


The network structure and parameters are built into the JSON format our associates at 
LANL are using and which our model is built to read. 


4.2 IEEE 13 Node Test Feeder Performance 

The IEEE RTS 13 node test feeder runs in a little over a second and returns an optimal 
solution objective function value of 0.00, a result our associates at LANL confirmed. 
Electrical behavior is captured in Table 4.3. 
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Table 4.3. IEEE 13 Node Test Feeder Edge Results 


Edge 

Type 

X ij 

Pij,k 

Qij,k 

634,633 

Transformer 

1 

1015 (a) 

381.91 (a) 




919 (b) 

515 (b) 




1132 (c) 

451 (c) 

633,632 

Line 

1 

1015 (a) 

381.91 (a) 




919 (b) 

515 (b) 




1132 (c) 

451 (c) 

632,650 

Line 

1 

0.00 

0.00 

632,645 

Line 

1 

400 (b) 

257 (b) 

645,646 

Line 

1 

230 (b) 

132 (b) 

632,dload 

Line 

1 

1015 (a) 

381.91 (a) 




519(b) 

258 (b) 




1132 (c) 

451 (c) 

dload,671 

Line 

1 

998 (a) 

371.91 (a) 




453 (b) 

220 (b) 




1015 (c) 

383 (c) 

671,684 

Line 

1 

128 (a) 

86 (a) 




170 (c) 

0.00 (c) 

684,611 

Line 

1 

170 (c) 

0.00 (c) 

684,652 

Line 

1 

128 (a) 

86 (a) 

671,680 

Line 

1 

0.00 

0.00 

671,692 

Line 

1 

485 (a) 

115.11 (a) 




68 (b) 

21.33 (b) 




460 (c) 

197 (c) 

692,675 

Line 

1 

485 (a) 

65.91 (a) 




68 (b) 

0.00 (b) 




290 (c) 

12(c) 


Our ACPF model demonstrates full functionality at this point, managing cycles, optimal 
switching, and returning expected operational results on a standardized test circuit. 
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Additionally, running the ACPF model against variations in our test networks (e.g., removing 
components, constraining electrical parameters, and varying loading amounts) produced 
electrical results we would expect for these variations. 


4.3 Expansion of IEEE 13 Node Test Feeder 

Our model is ultimately about performing analysis on distribution systems under the stress 
of disruption. With that in mind, we build an expanded version of the IEEE 13 node test 
feeder, essentially cloning the network, with “a” and “b” versions of each node and line. 
The two versions connect via a line from node 671a to 671b. This 3-phase line contains a 
switch which is open under normal operating conditions. This expanded IEEE 13 node test 
feeder is depicted in Figure 4.2 



The duplicated 13 node test feeder with “a” and “b” sides, connected by a line 
between nodes 671a and 671b. Electrical parameters for both sides are the same 
as the original IEEE 13 node test feeder, described in Tables 4.1 and 4.2. The 
cross connection line is 3-phase with Rjj and Xjj values copied from line 692-675. 
The Tij value is very low at 15 amps. 

Figure 4.2. Duplicated 13 Node Test Feeder 


When given to our ACPF model the duplicated IEEE 13 node test feeder sheds no load, 
an expected result since the single version met all load demands. In essence, we have 
parallel instances of the previous test network and each one saw commensurate electrical 
behavior, analogous to the results in Table 4.3 for both the “a” and “b” sides. Our next 
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step is to consider model performance and distribution system behavior in the presence of 
a disruption. 


Our chosen disruption is failure of the “b” side substation, essentially removing it from the 
network. Running the model with this disruption results in all “a” side load satisfied and the 
cross connection switch to close. Including the cross-connection, only a handful of lines 
on the “b” side have flow. Their values are captured in Table 4.4. 


Table 4.4. Duplicated IEEE 13 Node Test Feeder Edge Results 


Edge 

Type 

x ij 

Pij,k 

Qij,k 

67 la,67 lb 

Line 

1 

58.07 

58.07 

671b,dload_b 

Line 

1 

9.52 (a) 

38.82 (b) 

41.40(c) 

5.29 (a) 

13.90 (b) 

18.24 (c) 

Pij,k and %,k A 

lows the same for all phases unless exp] 

icitly listed. 


Due to relatively low thermal limit we imposed on the 671a - 671b connection, most of the 
“b” side goes without power and the model returns an objective function value of 4689.38. 
We are not surprised by this amount of load shed given the parameters we assigned the 
cross-connection. 


4.4 Model Performance Overview 

At this point the model has demonstrates the capability to analyze standardized test circuits 
in the presence of network disruptions. Electrical parameters behave as expected and 
operation of any distribution system can be optimized by the ACPF model in an effort to 
minimize load shed. We provide recommendations for follow-on work and additional model 
applications in Chapter 5. 
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CHAPTER 5: 

Conclusions and Follow-on Work 


5.1 Conclusions 

Working within a Python and Pyomo environment with CPLEX as a solver, our model 
demonstrates the ability to find an optimal network configuration and associated power flow 
for a 3-phase electrical distribution system. Our ACPF model can handle phase changes, 
varying and imbalanced loads, real and reactive demand, switch operation, and ensures the 
network is configured as a spanning tree to meet standard distribution system behavior. 

The novelty in our work is due to the focus on distribution systems. The model receives 
input from a JSON data structure used by the individuals at LANL working on infrastructure 
analysis. This allows for downstream use. Furthermore, the free and open-source nature of 
the Python programming language makes our work easily accessible for additional research. 

Alternating current is inherently non-linear. That said, the mathematics of the model are 
linearized, with the exception of constraint (4) which is non-linear. Testing and experi¬ 
mentation revealed that the non-linear structure of constraint (4) allowed for nonsensical 
solutions when implemented independently. As such, we identified the need to build and 
include constraint (4a). Doing so achieves correct results from the test networks and the 
IEEE RTS 13 node test feeder. 


5.2 Follow-on Work 

Our associates at LANL have a detailed and highly comprehensive dataset from a real 
northeastern utility provider. Our test networks and IEEE RTS 13 node test feeder are built 
in the same JSON format as the real-world data, facilitating model application to the real 
world data after some preliminary data cleaning. 

The real-world dataset contains 52 substations, 73,304 loads, and several hundred thousand 
nodes and arcs. There is interest in systematically triggering failure of substations and 
capturing network behavior after the ACPF model runs. This analysis would yield data for 
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subsequent vulnerability analysis as test and training sets. Furthermore, it would inform 
us regarding distribution system behavior in the presence of substation failure, a topic of 
interest to various government agencies. 

To facilitate future work investigating the behavior of distribution systems in the presence of 
adverse events, the model should be augmented to more easily disable part of the distribution 
system being analyzed. At present, the data itself needs to be modified to break a line or 
disable a substation. An update with the ability to pass in a component identification number 
as a function argument for disruption purposes would greatly simplify this process. 

Lastly, our ACPF model is built as a single Python function, allowing for easy incorporation 
into AD and DAD analysis. Analysis of this nature would yield insight to the system 
vulnerabilities of any distribution network and the corresponding hardening efforts of the 
greatest benefit. 
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